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Abstract 

We present an approximate method for calculating the electrostatic free energy of concentrated 
C^ protein solutions. Our method uses a cell model and accounts for both the coulomb energy and the 

entropic cost of Donnan salt partitioning. The former term is calculated by linearizing the Poisson- 
Boltzmann equation around a nonzero average potential, while the second term is calculated using 
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■^ a jellium approximation that is empirically modified to reproduce the dilute solution limit. When 
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combined with a short-ranged binding interaction, calculated using the mean spherical approxima- 



i-Q tion, our model reproduces osmotic pressure measurements of bovine serum albumin solutions. We 

Q also use our free energy to calculate the salt-dependent shift in the critical temperature of lysozyme 
• I— I 

"-^ solutions and show why the predicted salt partitioning between the dilute and dense phases has 

^^ proven experimentally elusive. 
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I. INTRODUCTION. 

Dense protein solutions are encountered in the cytoplasm and in vitro situations like 
pharmaceutical formulations, crystallization screens, and ultrafiltration [1-3]. In all these 
cases, the stability of the solution depends on sufficient electrostatic repulsion to overcome 
the short-range attraction of H-bond, hydrophobic, and van der Waals interactions. Elec- 
trostatic interactions are easily adjusted through changes in the pH or salt concentration 
providing a convenient experimental means to manipulate the phase behavior. However, 
they are difficult to model theoretically due to the long range nature of the coulomb force 
and the nonlinearity of salt screening. Because of this, most theoretical work on protein- 
protein interactions has focused on dilute solution properties where the electrostatic free 
energy is dominated by the coulomb energy and can be treated using two-body potentials 
[4-6]. We have recently shown that the electrostatic free energy of protein association is 
dominated by the change in salt ion entropy, which renders the net interaction strongly 
non-pairwise [7] . In this paper we extend these results in order to model electrostatic effects 
in non-ideal protein solutions. 

We test our theory against pH dependent measurements of the osmotic pressure and 
salt effects on the liquid-liquid phase separation of protein solutions. The osmotic pressure 
provides a direct test of the effective interparticle repulsion, while the latter phenomenon 
has attracted considerable attention due to the finding that fluctuations associated with 
the liquid-liquid critical point have been shown to accelerate crystal nucleation [8]. Liquid- 
liquid separation is analogous to the liquid-vapor coexistence in small molecules except that 
in proteins the liquid-liquid binodal is found entirely below the solid (crystal) solubility line 
due to the short range of the protein-protein attractive forces [9, 12]. Previous work has 
primarily focused on temperature as a means of controlling a phase behavior. However, 
temperature is a poor variable in protein systems because the accessible temperature range 
is limited to ~ 20% due to the freezing point of water and the thermal denaturation of the 
proteins. Our work shows how pH and salt can be used to manipulate the phase boundary 
into the accessible range. 



II. MODEL 

Our calculations are based on the following free energy of the ternary protein-salt-solvent 
system 

f{r]) = F{r])/N = f^v) + /att(r7) + fs.M + /couiomb(r7) (1) 

where f{ri) is the total free energy per protein as function of protein concentration r]. /he 
is the mixing entropy of spherical proteins interacting by a hard core potential, /att is the 
free energy due to short range attractions between the proteins, /gait is the mixing entropy 
of the salt ions, and /coulomb is the coulomb energy of the system. These terms are explicit 
functions oi rj = Na^/6V, the volume fraction occupied by the proteins, where a is the 
protein diameter and N/V is protein number density. We discuss each of these terms in 
detail below. 

A. Free energy of attractive spheres 

The terms /he and /att represent the free energy of a solution of attractive spheres. We 
adopt the attractive Yukawa potential to describe the short range interaction 

{oo r < a; 
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where 2 is a parameter describing the range of interaction, a is the hard sphere diameter, 
and e is the temperature dependent strength of interaction. Following previous work [13, 14], 
we take z = 4 reflecting the short range nature of the hydrophobic and H-bond interactions 
that dominate protein-protein attractions. 

Eq. 2 gives the binding energy of a two-body protein-protein interaction as a function 
of the center-to-center distance r. In order to derive macroscopic properties of the protein 
solution, we need to know the average binding energy per particle as a function of the protein 
concentration. We obtain this using the Mean Spherical Approximation [15, 16] the mean 
spherical approximation (MSA) is the method used to obtain the analytical solution of the 
radial distribution function of particles. The thermodynamic properties obtained from the 
MSA are in good agreement with results obtained from the Monte Carlo simulation. 



Within this approximation the binding energy density is [14, 17] 
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where X and Y are the variable defined in A9 and AlO. z is the parameter for range of 
interaction in the Yukawa potential and r] is the volume fraction of the proteins, ao, $0, 
J-'(X) and its first derivative are defined in Eqs. A2, A5, A12 and A13. This binding 
energy is partially offset by the loss of translational entropy at the high concentration. This 
contribution is given by the Carnahan-Starling expression [18] 
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hMlkBT = v'j^^r (4) 

Advantage of using MSA for our case is that, we get an expression for free energy density in 
terms of density or volume fraction 7] of protein. We can compute equation of state, osmotic 
pressure and many other thermodynamic quatities using this free energy density. MSA is 
based on the inverse temperature expansion of the free energy, hence it gives better result 
for higher temperatures. In our phase coexistence curve 6 and 7 we find our theoretical 
curves are more promising for large salt or when we have larger critical, temperautre. 

B. Electrostatic terms 

To solve for the electrostatic free energy, we adopt a cell model [19] in which each protein 
is surrounded by a spherical shell of solvent of radius h. The thickness of this solvent layer 
is chosen to reproduce the volume fraction occupied by the protein t] = {a/2b)^. We assume 
a protein solution in osmotic equilibrium with a reservoir of symmetric, monovalent salt 
of concentration Cg. The proteins carry a charge qCp, where Cp is the proton charge, that 
we take to be uniformly distributed over their spherical surface. The protein charge will 
perturb the salt ion concentration resulting in a local enrichment of counterions and a local 
depletion of colons. The extent of this ionic perturbation is a competition between the 
coulomb interaction of the salt ions with the electrostatic potential $ and the entropic cost 
of enriching/depleting the counterion/coion populations. These considerations are refiected 
in the electrostatic free energy /es = /coulomb + /salt where 

/coulomb = I y {V^fd^r (5) 

/sait/AJijT = / [c+ ln(c+/cs) -C+ + Cs] + [c_ ln(c_/cs) - c_ + Cs]d^r, (6) 
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where c± are the local cation/anion concentrations, V is the cell volume, V is the solvent 
accessible cell volume, and e is the local permeability which we take to be SOeo- The 
electrostatic potential can be expressed in terms of the charge distribution using the Poisson 
equation 

-V ■ (5V^(r)) = pp{r) + ep(c+(r) - c^{r)) (7) 

where pp is the charge density of the protein. The ion concentrations can be found by 
minimizing /es with respect to c± after integrating Eq. 5 by parts and applying Eq. 7. This 
resulting concentrations are 

c±(r) = c,e^^-*(")/^^^, (8) 

which, with Eq. 7, gives the well known Poisson-Boltzmann equation 

-V • {eV^ir)) = -ppir) - epC,(e-^-*W/'=^^ - e^vHr^/ksTy (g) 

Our strategy is to develop approximate solutions of Eq. 9 that can be used in Eqs. 5 
and 6 to calculate the free energy as a function of protein concentration, charge, and salt 
concentration. 

Within the cell geometry, the Poisson-Boltzmann equation has the approximate solution 

^(x) = (f){x) + 00 

^ c ( -'"''^-^' + f:!(i±A) _ t,„h^„ + ^„. (10) 

\ X X J 

where a = {cosh(j)Y^'^Ka/2, (3 = {cosh (pY^'^nb and C is a constant defined in the appendix. 
Eq. 10 is derived by linearizing the Poisson-Boltzmann equation around a reference po- 
tential 00 (00 = in the Debye-Huckel solution). The nonzero reference is necessary to 
model concentrated solutions of charged proteins, where the Donnan effect ensures that the 
potential never approaches zero. 

Eq. 10 can be used with Eq. 5 to determine the coulomb energy per particle. A 
straightforward integration yields (see appendix) 

/coulomb = -2C'(3' + —((/?- l)e-"+'^ + (/3 + l)e"-'^)2 

a 

-^((/3 + l)2e2("-/3) - (/3 - iYe~^i-~^^ - A{/3^ - l)a) (11) 

To compute /gait we first combine Eqs. 6 and 8 to obtain 
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This expression does not lend itself to direct integration, however, at high protein concentra- 
tions the potential varies weakly in the voids between proteins and it is reasonable to replace 
it with an average value (pix) — )■ (the jellium approximation) so that c± are constant [20]. 
Charge neutrality requires 

q = -Vion{c+ - c_) (13) 

= 2vionCsSmh{e(l)/kBT) (14) 

sinh"^ -— ^ — , (15) 
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where Vion, the solvent volume associated with each protein, is given by 

^'ion = Vpi 1 J , (16) 

where Vp = ira^/G is the volume of a single protein. 

Combining Eqs. 6 and 15 we find an expression for the salt entropy per protein at high 
protein concentrations 



fs.it/NknT = g(sinh-i ^ - v^l + ^^ + C) (17) 

where ^ = g/2fionCs. This expression is an excellent approximation for high concentration 
solutions, but fails in the dilute limit where it erroneously predicts that the entropic cost of 
the ion screening layers approaches zero. The problem can be traced to Eq. 16 which implies 
that the ion screening layer can become arbitrarily large. In reality, in dilute solutions the 
screening layer is confined to a shell with a thickness on the order of a Debye length. Because 
of this, the salt entropy will saturate at a minimum value when the concentration drops below 
a critical value rjQ. This behavior can be obtained from our model by a numerical integration 
of Eq. 6 with the potential given by 10. However, a more computationally efficient solution 
is to modify Eq. 16 to give the correct asymptotic behavior. This can be done with the 
following functional form 

%f = Vp It ^^T^ - 1 I • (18) 

Here n and rjQ are adjustable parameters that give optimal results when n = 5 and rjQ = 
(1 + 3.8//to")^ where /t^^ is the Debye screening length. Fig. 1 shows an excellent agreement 
in the salt entropy calculated using these two methods. Therefore, the remainder of the 
paper will utilize the more efficient effective water shell calculation of salt entropy. 
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FIG. 1: Comparison of two methods for calculating the salt entropy per particle. Solid lines show 
a numerical integration of Eq. 6 with the potential given by Eq. 10. Dashed lines show Eq. 17 
with the screening layer volume given by 18. cr/2 = 16nm, b = (o"/2 + 1.9/k), q = 10.3. The dashed 
lines have been shifted by an unimportant constant to facilitate comparison. 

III. RESULTS AND DISCUSSION 
A. Osmotic pressure 



As an initial test of our theory we compute the osmotic pressure of charged protein 
solutions 



U = 7] 



dl 
drj 



(19) 



Fig. 2 compares the theory to experimental measurements of bovine serum albumin [21]. We 
limit our comparison to pH values where the protein charge has been measured [22] , however 
the resulting charge values q = —20.2 (pH = 7.3), q = —9.1 (pH = 5.4), and q = 3.2 (pH 
= 4.6) cover a sufficiently large range to provide a rigorous test of the model. We use a 
protein radius of 2.4 nm [23] and the single fitting parameter is e = 1.22kBT = 2.98 kJ/mol. 
The agreement is generally good, with an average error of 6.5% for q = 3.2 and —9.1 and a 



larger 14% error (12% excluding the outlier at 100 g/1) at g = —20.2. 

In earlier work these data were modeled by fitting an effective protein volume for each 
pH [24]. These effective volumes ranged from slightly negative near the isoelectric point to 
four-fold greater than the actual volume when the protein carries a charge of ~ 60.0 [25]. 
Our modeling suggests that these trends can be explained by nonspecific protein binding 
competing with electrostatic monopole repulsion (mediated by the salt ions). 
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FIG. 2: Plot of osmotic pressure vs. protein concentration for three different protein charges, 
q = —20.2, q = —9.1 and q = 3.2 (solid lines). Points are data from [21] for three different pH 
values pH = 7.3, 5.4 and 4.6. The measured protein charges at these pH values are obtained from 
[22]. e = 2.98 kJ/mol is our fitting parameter. Temperature and salt concentration is same as in 
the experiment [21]. 

Fig. 3 shows predictions for the osmotic pressure as a function of salt for two values 
of the protein charge. Not surprisingly, greater salt concentrations result in lower osmotic 
pressures. This is explained by the fact that at higher salt concentrations the neutralizing 
counterions lie closer to the surface of the protein. Therefore, neighboring proteins will 
impose fewer constraints on the entropy of the screening layers. 
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FIG. 3: Plot of osmotic pressure vs. protein concentration for salt concentrations from 0.05M to 
l.OM. The two panels are for protein charges q = —20.2 (left) and q = —9.1 (right). 



Fig. 4 shows the component terms of the free energy as a function of the concentration. 
It is immediately apparent that the coulomb energy is minuscule in comparison with the 
other terms. This observation is somewhat misleading because the total free energy is a 
small residual from the addition of large terms with opposing signs. Therefore, the coulomb 
energy perturbation does give a significant quantitative improvement. On the other hand, if 
quantitative results are not required, a good estimate of pH and salt effects can be obtained 
from Eq. 17 while neglecting the more difficult to calculate coulomb energy. 

Another observation from Fig. 4 is that the repulsive terms are entropic in origin while 
the attractive terms are energetic. This has important consequences from the temperature 
dependence of the pressure. In Fig. 5 we plot the temperature dependence of the osmotic 
pressure for two systems of hard spheres. The first system interacts by the entropy dominated 
electrostatic interaction (Eqs. 11 and 17) while the second interacts by a repulsive Yukawa 
(Debye-Huckel) potential that is energy dominated (calculated using the MSA). We see that 
the temperature has qualitatively different effects on the two systems. In the Yukawa system 
the osmotic pressure actually decreases with temperature because the reduced Boltzmann 
weight given to the repulsive interaction at higher temperatures effectively frees volume for 
the spheres to explore. It is only with the entropy dominated electrostatic interaction that 
we recover the physically reasonable result that the osmotic pressure should increase with 
temperature [26]. 
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FIG. 4: Plot of the contributions to the free energy F/NKbT. The four curves are /he (solid), 
/coulomb (dashes), /att (dot-dash), /s^it (dots). The coulomb term is very small (inset) compared 
to other three terms. 

B. Liquid-liquid phase separation 

The solution free energy Eq. 1 contains a binding energy term that favors dense states 
competing with entropic terms that favor dilute states. This means that under some con- 
ditions the solution may phase separate in order to maximize the free energy contribution 
from these extreme states. The conditions for an equilibrium between two phases of unequal 
densities are 



/^, 



/i. 



/i. 



f^p 






(20) 
(21) 
(22) 



where /i' and /i" are the chemical potential of the dilute and dense phases, respectively, and 
the subscripts represent proteins, cations, and anions. The conditions for the salt ions are 
satisfied by Eq. 8 which captures both mixing entropy and coulomb energy contributions 
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FIG. 5: Plot of osmotic pressure (OP) 11 vs. protein concetration r/ for two different representations 
of the electrostatic interaction. Solid curves show the OP calculated for hard spheres with the 
entropy dominated interaction given by Eqs. 11 and 17. The three curves show the expected 
increase in OP with temperature. The dashed curves show the OP calculated for hard spheres 
with a Debye-Huckel type electrostatic interaction. This model shows a non-physical decrease in 
the OP with temperature. 

to the chemical potential. Eqs. 21 and 22 are, therefore, built into our free energy. Thus, 
the slope dF/dr] is an effective chemical potential for the protein that includes the effects 
of maintaining the salt equilibrium. 

The densities of coexisting phases are found by numerically searching for a line with two 
points tangent to F{rj) [27]. In Fig. 6 we plot the coexistence curve predicted by Eq. 1 
for three different NaCl concentrations (3%, 5% and 7% w/v) for a particle of diameter 
a = 3.2nm and charge q = 10. 3e. These parameters correspond to the pH 4.5 conditions 
used by Muschol and Rosenberger [28] to determine the liquid-liquid binodal curves for 
lysozyme. The single free parameter, e = 3.98 kJ/mol, has been set to reproduce the critical 
temperature under 7% salt. 
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It is immediately obvious that the theory fails to capture the width of the lysozyme 
coexistence curves. This is because proteins have much broader binodal curves than sys- 
tems of smooth spheres [11]. The reasons for this are not clear, but possible factors include 
anisotropic or directional binding [29-31], temperature dependence of the short range at- 
traction [13], asphericity [32], and the range of the attraction [14]. However, the theory 
does a reasonable job of reproducing the salt-dependent shift in the critical temperature, 
although as we observed with the osmotic pressure, the theory somewhat over-predicts the 
trend under the most extreme conditions (low salt or high charge). 
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FIG. 6: Plot of liquid-liquid coexistence curve for three different salt concentrations NaCl 3%, 
5% and 7%. Data points are experimental curve from [33]. We have good predictions for critical 
temperature. Other parameters protein charge q = 10.0, range of Yukawa interaction z = 4.0 and 
binding energy e = 3.98 kJ/mol. We kept binding energy as our fitting parameter and fixed all 
other parameters. 

A key feature of this theory is the partitioning of ions between the dilute and dense 
phases by the Donnan effect. However, experiments to observe this partitioning have yielded 
negative results [28]. Using Eq. 13 we can see why the partitioning was experimentally 
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FIG. 7: Plot of critical temperature Tc vs. salt concentration in Molar for protein charge q = 10.0 

elusive. At pH 4.5 lysozyme has a calculated charge of approximately Z = +10. If we take 
the protein volume to be 17 nm^ and assume that the dense phase has a density rj = 1/3, 
then for 0.5M salt conditions Eq. 15 gives a potential of eip/kBT = 0.46. From Eq. 8 we get 
a counterion concentration of 1.58 Cq within the solvent fraction of the dense phase. To get 
the apparent concentration within the total phase we multiply by 1 — r^ to account for the 
excluded volume of the protein. We find that the apparent counterion concentration in the 
dense phase is nearly unchanged from the bulk value. However, this cancelation is entirely 
a coincidence of the conditions used in Ref [28] . A similar calculation for the coions yields 
an apparent concentration of 0.4co, however, this may have been missed due to difficulties 
associated with the sodium assay [M. Muschol personal communication]. 

IV. CONCLUSION 

We have presented an approximate method for calculating electrostatic effects in dense 
protein solutions. Our method captures the non-pairwise nature of electrostatic interactions 
at high concentrations, yet requires only a single particle calculation within a cell model. 



13 



Such methods may enable the rational manipulation of protein solution behavior and phase 
diagrams using pH and salt concentration as control parameters. 
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Appendix A: Mean Spherical Approximation MSA 

The attractive term fatt in Eq. 1 can be given by the MSA [15] for a Yukawa fluid. As a 
second-order type theory, the MSA yield 



^ oo 
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By doing some algebric manipulation the free energy takes the simple form 
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Appendix B: Linearized Potential 

Here we solve for the ion distributions using the cell approximation. We model the protein 
as a charged sphere embedded in a spherical cavity of solvent with radius h = ar]^^^^. To do 
this we need to solve the Poisson-Boltzmann equation 



eVl^{f) = -pf{f) - eiV^(coe-^*(^^/'=-^ - coe'"^"^/'-^). 



(Bl) 



for the electrostatic potential \E'. Here pf is the fixed charge distribution on the protein, cq 
is the bulk salt concentration, A^^ is Avogadro's constant, and e ~ 80eo is permeability of 
water. First, we put the PB equation in dimensionless form 

Vj* = sinh($) (B2) 

where $ = Cp^ /ksT is the dimensionless potential, y = /tr is a dimensionless length, and 

.^ ^ ^. (B3) 

Now we linearize the potential around the local potential $(|/) = 0(|/) + 0o, so we have 



Vy<l> = sinh(0 + 0o) 



) cosh 00 + sinh 0o 
) + tanh0o, 



where the new length variable is 
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The solution to Eq. B6 is 



(x) = A h B tanh 0o, (B9) 



where the constants A and B are determined by the boundary conditions. The boundary 
conditions are that at the inner sphere r = a the electric field is equal to that of the bare 
macroion, and at the outer sphere r = b the electric field vanishes (due to charge neutrality). 
After scaling the sphere radii 

a = a/ cosh0o'^fl (BIO) 

/3= ^/cosh(f)oKb (Bll) 

the inner boundary condition becomes 

= -A— (l + -)+ B- (l--), (B13) 

a \ a) a \ a) 

where the dimensionless electric field is 

gcp/ty/cosh^o 
^° = ^-aek^Ta^ ' ^^^^^ 

where qe^ is the charge on the central sphere. At the outer sphere boundary we have 



from which we derive 



i.^.4e-»l±i^. ,B16) 



Combining Eqs. B13 and B16 we find 

A[-e-"(a + l)(/3 - 1) + e"-2^(a - l)(/3 + 1)] = -a'Eo(/3 - !)• (B17) 

So the two constants are 

. ^ ^'E,ef^{ii-\) 

e-"+/3(a + l)(/3-l)-e"-/3(a-l)(/3 + l) ^ ^ 

B = -^E,e-^iP + l) 

e-"+/3(a + l)(/3-l)-e"-/3(a-l)(/3 + l)' ^ ' 
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The total potential is then 
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so the constant used in the text is 
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The coulomb energy is obtained as follows 
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where Ai = C{(3 - 1) exp(/3) and 5i = C(/3 + 1) exp(-/3). 
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